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Abstract. Aims: morpholgies, number and energy distributions of Cosmological Shock 
Waves from a set of ENZO cosmological simulations are produced, along with a study of 
the connection with Cosmic Rays processes in different environments. Method: we perform 
cosmological simulations with the public release of the PPM code ENZO, adopt a simple 
and physically motivated numerical setup to follow the evolution of cosmic structures at the 
resolution of \25kpc per cell, and characterise shocks with a new post processing scheme. 
Results: we estimate the efficency of the acceleration of Cosmic Ray particles and present 
the first comparison of our results with existing limits from observations of galaxy clusters. 
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1. Introduction 

Galaxy clusters store up to several 10 63 ergs 
in the form of hot baryonic matter, due to 
the action of shock-heatin g processes in larg e 
scale structures formation dZel'Dovicri ll970). 
Detecting shocks in Large Scale Structures 
(LSS) is still observationally challenging since 
they should more frequently develop in pe- 
ripheral, low X emitting regions of clusters, 
due to the drop in the sound speed there (e.g. 
iMarkevitch & Vikhhnin1l2007l) . Shocks are im- 
portant not only to understand the heating of 
the ICM but also because they are expected to 
be efficien t accelerators of supra-thermal parti- 
cles (e.g. lSarazinlll999l) . which can then be ad- 
vect ed and accumulated inside galaxy cluster s 
(e.g. IVolk et al.lfl996l iBerezinskv etai1ll997l) . 
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Non thermal activity in galaxy clusters related 
to the presence of Cosmic Ray (CR) electrons 
and ~ yi/G magnetic field is proved by radio 
observations, which show synchrotron emis- 
sion in a fra ction of merging clusters (e.g. 
iFerettil 120051) . in form of Radio Haloes (at 
the cluster center) and Radio Relics (elon- 
gated and at the cluster periphery). CR pro- 
tons are expected to give y-Ray emissions 
from galaxy clusters via the decay of tt° gen- 
erated during proton-proton collisions in the 
intra galactic medium. Still, only upper lim- 
its of this y-Ra y emissions have b een ob- 
tained so far (e.g. iReimer et al.ll2003l) . Recent 
numerical works claimed that an efficient CR 
protons acceleration can oc cur in larg e scale 
shocks dMiniati et all l200 ll iRyu et all 120031: 
iPfrommer et all 120061) . However, the identifi- 
cation and characterisation of shocks, as well 
as the calculation of the energy injected in the 
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form of CR, remain challenging due to the un- 
certainties in the numerical schemes and due to 
our ignorance of the efficency of the accelera- 
tion of CRs at relatively weak shocks. 

In this paper we follow th e appr oach of 
the seminal paper bv lRvu et alj d2003l) . study- 
ing the shock wave patterns in LSS. Shocks 
are characterised in a post-processing phase 
with a new scheme, which evaluates shocks 
Mach number by analysing the jumps in the 
velocity varaible. Estimates of the level of 
CR injections at these shocks are provided 
and compared to present day observational 
upper limits. More detailed discussions and 
presentations of the results can be found in 
Vazza, Brunetti & Gheller (2008 submitted to 
MNRAS, hereafter VBG). 

2. Numerical Code - ENZO. 

The simulations described in this paper were 
perfo rmed with the cosm ological code ENZO 
(e.g. lO'Shea et all |2004 . ENZO couples an 
N-body particle-mesh solver for Dark Matter 
with an adaptive mesh method for ideal flui- 
dynamics. The hydrodynamical solver is based 
on the the Piecewise Para bolic Method (PPM, 
IColella & Woodwardll 19841) . which is a higher 
order extension of Godunov's shock captur- 
ing method. It is at least second-order accu- 
rate in space (up to the fourth-order, in the 
case of smooth flows and small timesteps) and 
second-order accurate in time, and it is thus a 
highly suitable hydro method to study shock 
pattern. 

3. Cosmological Simulations. 

We adopt the "Concordance" model, with den- 
sity parameters Qo = 1.0, Q.bm = 0.044, 
£Idm - 0.226, D.a = 0.73, Hubble parame- 
ter h = 0.71, a power spectrum produ ced ac- 
cording to the lEisenstein & Hul (1999) fitting 
formulas with a primordial spectrum normal- 
ization <t 8 = 0.94, and an initial redshift of 
z = 50. In order to have a large number stas- 
tics of massive galaxy clusters, we collect sev- 
eral boxes in order to produce a final equivalent 
volume of * (IQQMpc/h)^ at the fixed numeri- 
cal resolution of 125 kpc, Our fiducial model 



here is an ensamble of non-radiative simula- 
tions with a a post-processing treatment of cos- 
mic reionizaation, desi gned to reproduce the 
lHaardt & Madaul dl996h model. 



4. The Velocity-Jump Method 

The crossing of a shock in a simulated volume 
leaves its imprint as a jump in all the thermo- 
dynamical variables, which can be inverted to 
evaluate the shock Mach number, M, by means 
of the standard Rankine-Hugoniot jump con- 
ditions. The relationship between the jumps in 
the velocity field, Av, and M for an idealized 
shock wave running in an imperturbed medium 
follows from the conservation of momentum 
and density across the shock, and transforming 
the velocities from the shock frame to the Lab 
frame we get: 



3 1 - M 2 
Av — — v s - — . 

4 M 2 



(1) 



where v s = Mc s and c s is the sound veloc- 
ity computed in pre-shocked cells. 

The procedure adopted is the following: 
1) we consider only "candidate" shocked cells 
with a negative 3-D velocity divergence; 2) 
in the case of two adiacent candidate shocked 
cells, the one with the minimum 3-D diver- 
gence is considered as the post-shock region; 

3) we perfrom 1-D scan along each axes mea- 
suring all Av A V Z jumps across neighbours cells; 

4) we measure the shock Mach number along 
each coordinate according to Eq[T] where v s 
is calculated from the temperature of the pre- 
shock cell; 5) the total Mach number for the 
shocked cell is finally calculated as M - (M 2 + 
M 2 + M 2 ) 1 ' 2 . 

As in other methods relying on a post- 
processing of the simulated output, this method 
has its major source of uncertainty in the as- 
sumption of an unperturbed velocity field prior 
to the passage of the shock. However, in the 
case of cosmological simulations the situa- 
tion is more complex due to the chaotic pat- 
tern of velocity and temperature fluctuations 
which develops d uring LSS formation (e.g. 
iDolag etal.1 120051) . The uncertainties deriving 
from that are discussed in VBG, where we 
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Fig. 1. 2-dimensional slice of the recon- 
structed Mach numbers of shocks for a box of 
side SOMpc. 

claim that in the case of dense cosmological 
regions the efficency of the Velocity Jump (VJ) 
method in characterising shocks is similar to 
that of methods based on temperature jumps 
across cells within clusters, and becomes bet- 
ter in lower density regions (i.e. filaments and 
voids). 

5. RESULTS 

In this Section we present the main results ob- 
tained for shocks in LSS at z = with the 
VJ method. Additional results concerning the 
time evolution of shocks in the cosmological 
volume can be found in VBG. 

5.1. Maps and Number distributions. 

Fig Q] shows the distribution of detected shocks 
in a I25kpc cut of a region of side SOMpc, 
taken at z = 0. Roughly a ~ 15 per cent of 
the simulated volume hosts shocks at present 
epoch, with the percentage of shocked cells 
increasing in denser environments. Overall, 
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Fig. 2. Distribution of shocks Mach number for 
the whole fiducial data-set (black line) and dis- 
tributions for subsamples in the cosmic envi- 
ronment (colour lines). 

th e picture is pretty s imilar to early results 
bv lMiniati et al.l d2000h . whith filamentary and 
sheet-like shocks developing at the interfaces 
of cosmic filaments and sheets, and with 
shocks surrounding galaxy clusters showing 
spherically shaped boundaries at a typical dis- 
tance of about 1 - 2R V i r from clusters cen- 
ter. Internal merger shocks are more irregu- 
lar and weak, M < 2, while slightly stronger 
shocks are only episodically found within clus- 
ters in case of merger events. An issue which 
is still poorly addressed in the literature is the 
Mach number distribution of shocks in numer- 
ical simulations. Figj2] shows the distribution 
of shocked cells found with the VJ method at 
z — 0, for the total volume and in different cos- 
mic environments. All distributions are steep, 
with a ~ -3.5 (a = d log N(M)/d log M) for 
the whole volume, and a < -6 inside the 
virial radius of galaxy clusters. The majority of 
shocks are weak, with their distribution show- 
ing everywhere a peak at M ~ 1.5 and a mono- 
tonic decrease at larger M. 

5.2. Thermal Energy Flux in Shocks. 

A shock wave thermalises the post-shock re- 
gion according to Rankine-Hugoniot jumps 
conditions, which relate the flux of the kinetic 
energy crossing the shock and the resulting 
thermal flux in the post-shock region, /,/,. We 
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Fig. 3. Left:3-D rendering of the dissipated energy flux for a region of 80Mpc per side. Color 
coding goes from blue (f,h ~ 10 i3 erg/s) to yellow (f,i, ~ 10 3& erg/s) to red (/,/, > 10 41 erg /s). 
Right: energy ratio between injected CR energy flux and thermal energy flux in shock waves, for 
a slice crossing the center of the same two clusters in left panel. 
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Fig. 4. Distribution of the thermalised energy 
flux at different overdensity bins, for the whole 
AD 125 and normalized to a comoving volume 
of (IMpc/h) 3 . 



follow the formal ism of iRyu et al.1 (l2003h and 
iKang et all d2007l) to calculate f,k at a given M. 

Left panel in Figj3] shows a 3-D render- 
ing of the thermal energy flux through shocks 



Q, while Figj4] depicts the distribution of f,i, 
at shocked cells, for different cosmic environ- 
ments. 

About 70 per cent of the total thermal en- 
ergy flux from shocks comes from the virial re- 
gion of galaxy clusters and the bulk of the ther- 
malisation is done by shocks with M « 1.8-2. 
These relatively weak shocks are also responsi- 
ble for the thermalisation in lower density envi- 
ronments, but in these regions a sizeable frac- 
tion of the thermal energy flux is injected at 
larger Mach numbers. 

the total thermalised 
in our simulations is 
finding, our distribu- 
find a, i, « -3 (with 
M a ' h ), which should 



We find that although 
energy per cubic Mpc/h 
consistent with previous 
tion of /,/, is steeper: we 
a t h taken as f t h(M)M oc 
be co mpared to a,/, « -1 
d2003|) and to a th x -2 



.5 to -2 in Rvu et al 
5 in Pfrommer et al 



1 The 3-D rendering is g enerated with of the 
visualization tool VISIVO jBecciani et"al] 120071 . 
http://visivo.cineca.it) 
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5.3. Cosmic Rays acceleration. 

The injection and acceleration of Cosmic Rays 
at shocks is a complex process, where sev- 
eral still unknown quantities play a major 
role (e.g. iBlasil 120041 for a review). In the 
case of a low level of energy in form of 
CR, it is customary to describe the accel- 
eration according to the diff usive shock ac- 



celeration (PSA) theory (e. g. Drurv & Voelkl 



11981b iBlandford & OstrikerlTl978h . This the 
ory applies when particles can be described 
by a diffusion-convection equation across the 
shock, but it fails in case of strong shocks, 
where the pressure in the acceler ated CR back 
reacts on the shock itself (e.g., lEUison et alJ 
fl99l . 

In order to have a straightforward compar- 
ison with other numerical papers on the is- 
sue, we estimate the ratio between the energy 
flux trough a shock and the energy flux which 
is channelled into CR acceleration, feu, by 
means of a simple parameter, t](M) = fcit/fy 
dKang et al.ll2007l for an analytical expression 
of T](M)), which depends on the Mach num- 
ber only. Fig|3](it right panel) shows the spa- 
tial distribution of fcit/fa at shocks for a cut 
taken across the same region of the left panel. 
The highest values of fcR I J "th are found at the 
interface layers of filaments or in the outer- 
most regions of galaxy clusters, where a sub- 
stantial population of relatively strong shocks 
is present. On the other hand the lower val- 
ues are typically found inside galaxy clusters, 
where the Mach number distribution is steep 
and strong shocks are extremely rare; here 
fcu/fih ^ 0.10 is measured. 

As in the case of the thermal energy, we 
find that the distribution of the energy in CR 
at shocks is steeper than that reported in other 
works, with the bulk of the CR injection tak- 
ing place at M ~ 2 - 2.5, at all cosmic envi- 
ronments. Since_weusean approach equivalent 
to that in lRvu et akl d2003l) to evaluate the CR 
acceleration, this difference is likely related to 
our different shock detecting scheme, and ex- 
pecially to our modeling of the re-ionization 
process in the simulations (see VBG for a de- 
tailed discussion). 



A compar i son with the results in 
iPfrommer et alJ (120061) is more difficult since 
these authors use a Lagrangian Smoothed 
Particles Hydrodynamics code which also 
include a different approach in the calculation 
of CR dynamics. The overall distribution of 
the energy flux inje cted in CR reported in 
IPfrommer et al.l (120061) has a slope cicr ~ —1.8 
which has to be compared with the value of 
<xcr ~ _ 2.2 that we find in our simulations. 

5.4. Shocks in Galaxy Clusters. 

In this Section we studied shocks in four rep- 
resentative galaxy clusters with masses above 
5 • 1O 14 M and with different dynamical states: 
a cluster in a relaxed state, a system with 
an ongoin merger with a smaller subclump, 
a system approaching a major merger, and a 
post-merger system (2 Gyr after the close en- 
counter). 

Fig|5]reports the distribution of fa with M 
in shocked cells within 1 /?,,,> from the clus- 
ter centers; the distributions in the four clus- 
ters were normalized to the volume of the most 
massive system (a sphere of radius ~ 3>Mpc). 
A very steep distribution of the energy flux 
trough shocks is found in the relaxed and in 
the minor merger case, while in the case of the 
ongoing merger and in the post merger case 
the distributions also show tail of higher Mach 
numbers. Inside /? V! > no shocks with M > 3 
are detected, except for a few in the case of 
the post-merger system, and this is in line 
with X-r ay observations of real merging clus- 
ters (e.g. iMarkevitch & Vikhlininl 120071) . We 
remark that our findings are in line with ex- 
pectations from semi-analytical treatment of 
shocks in virialzied merging galaxy clusters 
dGabici&Blasill2003h . On the other hand, we 
find significantly less strong shocks than in 
Pfrom mer et al.l (120071) . In this case the differ- 
ences are due to the different numerical scheme 
and to the procedure adopted to characterise 
shocks, and highlight the level of present ucer- 
tainty in this issue. 

The radial behaviour of the fculfa for 
these four clust ers is reported in Figj6] where 
we adopted the lKang et al.l(l2007l) formulation. 
Inside the virial radius we do not find any rele- 
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Fig. 5. Distribution of thermalised fluxes for 
the four different galaxy clusters presented in 
the text. 
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Fig. 6. Volume averaged profiles of the injec- 
tion efficency, fcit/fth> f° r me f° ur galaxy clus- 
ters. 



vant difference among our clusters, and an av- 
erage value of fcit/fth < 10 is found. This 
is because shocks crossing the innermost re- 
gions are weak, independently of the dynam- 
ical status of the clusters. Although we do not 
follow the advection of CR with cosmic time, 
and this makes a comparison with observations 
more challenging, our estimates appear in line 
with existing upper limits from Radio observa- 
tions for uG magneti c field in galaxy clusters 
dBrunetti et all 120071) . as in our case the spec- 
trum of CR is steeper than previous work (with 
an average slope of 5 ~ -3.5, see VBG for a 
detailed discussion). 



6. Discussion and Conclusion. 

In this paper we have investigated the proper- 
ties of shock waves in LSS simulations pro- 
duced with the cosmological code ENZO and 
by means of a new detection scheme to mea- 
sure the shocks, based on velocity jumps. 

In the following we summarize the main re- 
sults obtained: 

- we detect morphologies of LSS shock- 
patterns which are qualitatively in agree- 
ment with previou s numerical works (e.g. 
MiniatTetdj2001), with strong shocks en- 
veloping filaments and sheets of matter, 
and weaker shocks hosted inside galaxy 
clusters. 

- We measure the number distribution of 
shocks as a function of the Mach in all cos- 
mological environments. The bulk of cos- 
mological shocks is made by weak M < 2 
shocks and their distribution can be grossly 
described by a steep power law N(M) oc 
M a , with a « -3.5. In the case of galaxy 
clusters a « -6, demonstrating the increas- 
ing rarity of strong shocks in these denser 
(and hotter) r egions. 



- Following Ryu et all (120031) and 



Kang_etaLl 12007) we calculate the 
energy rate dissipated in thermal energy at 
shocks, finding that roughly a 70 percent 
of the thermal dissipation in the whole 
volume happens inside galaxy clusters, 
at an average Mach number of M « 2. 
Although in qualitative agreement with 
previous studies, the energy distributions 
we measure in all enviro nments are s teeper 
than t hose obtained by iRvu et al.1 (120031) 
and by lPfrommer et all (120061) . 

- We calculate the efficeny of CR accel- 
eration for our simulations. Also in this 
case our results are in qualitatively in 
line with previous findings, although our 
energy distributions with M ar e sistem- 
aticalj y steeper than those in iRyu et all 
(2003) and slightly s teeper than those in 
iPfrommer et all d2006l) . 

- We report on the properties of shocked 
cells propagating in four representative 
galaxy clusters of our sample. The aver- 
age Mach number of shocks within li? VJ > 
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is M m 1.5 and this is in line with semi- 
analytical studies dealing with mergers o f 
virialised systems (Gabi ci & Blasil 120031) . 
Also, the rariry of stronger shocks (M > 
2 - 3) is found in line with the the rarity of 
shocks det ected so far by X-ray observ a- 
tions fe.g. lMarkevitch & Vikhlininll2007h . 

F. V. thanks K. Dolag, D. Ryu, H. Kang, 
G. Tormen, L. Moscardini, S. Giacintucci and 
R. Brunino for useful discussions and helps. 
We acknowledge partial support through grant 
ASI-INAF 1/088/06/0, and the usage of com- 
putational resources under the CINECA-INAF 
agreement. 
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